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Abstract 

A  two-dimensional,  non-isothermal,  two-phase,  multicomponent  model  is  presented  for  the  cathode  of  a  PEM  fuel  cell,  which  can  be  applied  using 
both  conventional  and  interdigitated  gas  distributors.  Gas  and  liquid  transport  are  described  by  an  extended  Darcy’s  law  using  a  continuum  approach. 
The  most  important  material  properties  like  wettability,  permeability,  porosity  and  tortuosity  are  determined  experimentally.  For  the  determination 
of  the  capillary  pressure-saturation  data,  a  method  based  on  mercury  intrusion  porosimetry  has  been  applied.  The  relative  permeability  was  then 
derived  from  this  data.  It  turned  out  however,  that  a  more  detailed  description  of  the  partial  wetting  behaviour  of  hydrophobic  electrodes  is  of 
decisive  importance.  The  two-dimensional  distribution  of  the  reactants  and  products  is  obtained  and  depends  strongly  on  the  flow  configuration  and 
the  operational  conditions.  The  model  results  were  compared  to  experimental  data  and  qualitative  information  on  the  effects  of  various  operating 
conditions  and  on  material  properties  were  obtained. 

©  2006  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Polymer  electrolyte  membrane  (PEM)  fuel  cells  with  hydro¬ 
gen  or  methanol  as  fuel  are  very  promising  for  applications  in 
electrical  vehicles  and  portable  power  sources  and  have  there¬ 
fore  been  intensively  investigated  during  the  past  decade.  Their 
low  operating  temperature  in  the  range  from  ambient  tempera¬ 
ture  up  to  approximately  100  °C  is  one  of  their  most  important 
advantages  for  mobile  and  portable  applications  in  comparison 
with  other  fuel  cell  types.  However,  further  efforts  have  to  be 
made  in  order  to  improve  the  performance,  reduce  the  costs,  and 
thus  make  these  systems  competitive  with  conventional  energy 
supplying  systems. 


Abbreviations:  CCD,  charge  coupled  device;  DS,  double  sided;  GDL,  gas 
diffusion  layer;  HOR,  hydrogen  oxidation  reaction;  ME  A,  membrane-electrode- 
assembly;  ORR,  oxygen  reduction  reaction;  PEM,  polymer  electrolyte  mem¬ 
brane;  SEM,  scanning  electron  microscopy 
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The  PEM  fuel  cell  consists  of  a  proton  conducting  poly¬ 
mer  membrane  situated  between  two  porous  catalytic  electrodes. 
Each  electrode  consists  of  a  porous  gas  diffusion  layer  (GDL) 
or  backing,  and  a  reactive  layer  containing  a  platinum  or  plat¬ 
inum  alloy  catalyst.  This  membrane-electrode-assembly  (MEA) 
is  then  placed  between  two  cell  plates  which  provide  electrical 
connection  and  have  flow  fields  that  direct  the  reactants  on  the 
surface  of  the  reactive  layer  and  also  permit  the  removal  of  the 
products  from  the  cell. 

The  chemical  reactions  in  a  hydrogen  fuel  cell  are: 

•  Anode: 


H2  -»  2H+  +  2e“ 

(1) 

Cathode: 

±02  +  2H+  +  2e“  -*  H20 

(2) 

Overall: 

h2  +  ^02  H20 

(3) 
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Nomenclature 

A  i  empirical  constant  in  Eq.  (27) 

Aj  empirical  constant  in  Eq.  (28) 

B\  empirical  constant  in  Eq.  (27) 

B2  empirical  constant  in  Eq.  (28) 

cs  heat  capacity  of  the  solid  matrix  (Jkg-1  K-1) 

Ci  empirical  constant  in  Eq.  (27) 

C2  empirical  constant  in  Eq.  (28) 

Z)eff  effective  diffusion  coefficient  (m2  s-  1 ) 

Djj  binary  diffusion  coefficient  (nr  s-1) 

Dpm  effective  diffusion  coefficient  in  porous  medium 
(m2  s-1) 

D\  empirical  constant  in  Eq.  (27) 

D2  empirical  constant  in  Eq.  (28) 

£2  empirical  constant  in  Eq.  (28) 

/v  catalytic  surface  increasing  factor 

F  Faraday’s  constant,  96484.56  C  mol- 1 
g  gravity  vector  (m  s-2) 

A Gq  free  activation  enthalpy  (J  mol-1) 
h  enthalpy  (kJ  kg- 1 ) 

i  current  density  (A  cm-2) 

ir(f  exchange  current  density  ORR  (A  cm-2) 

J  diffusive  molar  flow  (mol  s-1  m-2) 

kT  relative  permeability 

K  absolute  permeability  (m2) 

N  molar  flux  (mol  m-2  s- 1 ) 

p  pressure  (Pa) 

pc  capillary  pressure  (Pa) 

q  source  or  sink  term  (mols-1m-3),  heat  flux 

(J  s-1  m-2) 

qh  energy  source  or  sink  term  (J  s-1  m-3) 

r  radius  (m) 

R  universal  gas  constant  (8.3 14  J  mol- 1  K- 1 ) 

Rspec  specific  resistance  (£2  m2) 

S  phase  saturation 

Sa  phase  saturation 

t  transport  number,  time 

T  temperature  (K) 

u  internal  energy  (kJ  kg- 1 ) 

f7Ceii  cell  voltage  (V) 

Uoc  open  circuit  voltage  (V) 

Ulev  reversible  cell  voltage  (V) 

Uth  thermoneutral  cell  voltage  (V) 

v  Darcy  velocity  (ms-1) 

x  molar  fraction 

Z  number  of  transferred  electrons 

Greek  symbols 

a  fluid  phase,  load  transfer  number  in  Eq.  (12) 

<5  thickness  (m) 

q  overvoltage  (V) 

bequ  equilibrium  overvoltage  (V) 

Ptrans  transfer  overvoltage  (V) 

9  contact  angle 


X  thermal  conductivity  (W  K-1  m-1) 

/x  dynamic  viscosity  (kg  m- 1  s- 1 ) 

Vi  stoichiometric  coefficient 

pa  phase  density  (kg  m3) 

a  surface  tension  (N  m- 1 ) 

r  tortuosity 

4>  porosity 

Subscripts  and  superscripts 
a  anode 

c  cathode 

g  gas  phase 

HtO  water  component 

inlet  cell  inlet  conditions 

K  component 

M  membrane 

N2  nitrogen  component 

O2  oxygen  component 

pm  porous  medium 

r  residual 

R  reaction 

ref  reference  conditions  ORR 

w  liquid  water  phase 


Potential  losses  under  electrical  current  occur  due  to  kinetic  ac¬ 
tivation  resistances,  ohmic  losses  and  mass  transfer  limitations, 
as  can  be  recognized  from  the  form  of  a  typical  polarization 
curve  [1].  It  is  an  accepted  fact  that  the  cathode  of  the  H2- 
PEM  fuel  cell  is  the  performance-limiting  component  due  to 
the  slower  kinetics  of  the  oxygen  reduction  reaction  [2]  and  the 
mass-transfer  limitations  for  oxygen.  They  are  caused  by  the 
stagnant  nitrogen  (if  operated  with  air)  and  the  counter  diffu¬ 
sion  of  the  water  vapor  formed  and  are  enhanced  by  the  liquid 
water  generated  by  electrochemical  reaction  and  electro-osmotic 
drag,  which  accumulates  in  the  porous  catalyst  and  gas  diffusion 
layers. 

The  water  management  in  a  PEM  fuel  cell  is  difficult  due  to 
the  fact  that  the  membrane  has  to  be  humidified  to  be  able  to 
transport  the  protons  from  anode  to  cathode.  On  the  other  side, 
water  accumulating  in  the  cathode  electrode  hinders  the  oxygen 
transport  to  the  catalyst  layer. 

The  reactive  gases,  typically  hydrogen  and  air,  normally  have 
to  be  humidified  to  avoid  dry-out  of  the  membrane.  Thus,  a 
gas-liquid,  two-phase  flow  is  occurring  within  the  flow  field 
at  the  cathode,  particularly  at  high  current  densities  when  the 
amount  of  water  from  reaction  and  electro-osmosis  are  higher. 
The  water  phase  may  flood  part  of  the  electrode  pores  or  the 
channel,  hence  blocking  the  oxygen  supply  to  the  reactive  layer. 
Liquid  water  also  covers  the  catalyst  sites,  thus  inhibiting  the 
oxygen  reduction  reaction  (ORR). 

The  interaction  of  material  properties,  cell  and  electrode  de¬ 
sign,  and  the  operating  conditions  affect  the  water  management 
and  thus  the  performance  of  the  cell.  In  practice  the  backing 
is  made  of  conducting  carbon  material,  which  is  treated  with 
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Teflon  to  render  it  hydrophobic  and  enhance  the  liquid  water  re¬ 
moval.  Moreover,  new  types  of  flow  fields  have  been  developed 
to  improve  the  gas  flow  and  water  transport  [3]. 

Due  to  the  small  geometry,  the  actual  reactant  and  product 
concentration  profiles  in  the  catalyst  and  gas  diffusion  layers 
are  very  difficult  to  measure.  This  has  prompted  researchers 
to  develop  mathematical  models  to  provide  qualitative  insight 
into  the  processes  involved  in  the  fuel  cell.  The  modeling  has 
been  specially  focused  on  the  cathode,  due  to  its  importance  as 
performance-limiting  component. 

Previous  studies  were  mostly  based  on  a  simple  (gas)  phase 
assumption,  with  the  presence  of  liquid  water  being  neglected. 
Many  of  these  models  are  one-dimensional  (dimension  normal 
to  the  reactive  catalyst  surface),  e.g.,  [2,4],  and  also  in  1  +  1 
or  two  dimensions,  e.g.,  [5-10].  Recently,  several  studies  on 
two-phase,  multidimensional  transport  in  PEMFC  have  been 
made  [11-20].  Only  few  investigations  [12,17,18]  include  the 
energy  equation  to  predict  the  effect  of  heat  produced  by  the 
electrochemical  reactions  on  fuel  cell  performance.  Shimpalee 
et  al.  [12],  for  example,  treat  the  liquid  water  as  a  component 
of  the  gas  mixture  and  not  as  a  separated  phase.  Moreover, 
all  these  models  assume  hydrophilic  properties  of  the  gas 
diffusion  layers.  [20,21]  used  a  Leverett  approach  and  extended 
it  to  obtain  hydrophobic  data  for  use  in  their  one-dimensional 
models. 

For  the  typically  used  hydrophobic  carbon  gas  diffusion  ma¬ 
terials  there  is  a  lack  of  property  data.  Most  recently,  [22]  char¬ 
acterized  some  commercial  gas  diffusion  materials,  reporting 
data  like  contact  angle,  permeability,  porosity  and  conductivity. 
No  measured  capillary  pressures  have  been  published  for  those 
materials,  mainly  due  to  the  difficulty  in  applying  conventional 
methods,  e.g.,  see  [23],  to  these  small  dimensioned  materials. 

In  the  present  study  we  propose  a  simple  method  to  determine 
the  capillary  pressure-saturation  relationship  for  the  usual  car¬ 
bon  diffusion  layer  materials  and  implement  these  data  in  our 
two-dimensional,  non-isothermal,  two-phase  model.  We  com¬ 
pare  the  results  with  measured  polarization  curves  and  use  our 
model  to  predict  the  influence  of  some  parameters,  like  inlet 
stream  humidity  and  capillarity,  on  the  liquid  water  distribution 
and  cell  performance. 

In  the  following,  the  governing  equations  and  boundary  con¬ 
ditions  used  in  the  model  are  presented  in  Section  2.  The  applied 
numerical  discretization  technique,  given  in  Section  3,  follows. 
Section  4  shows  the  experimentally  determined  modeling  pa¬ 
rameters  followed  by  the  results  and  discussions  in  Section  5 
and  the  summary  and  conclusions  of  the  work,  which  are  given 
in  Section  6. 

2.  Mathematical  model 

2.1.  Model  description  and  assumptions 

Fig.  1  shows  a  schematic  representation  of  the  domain  mod¬ 
eled.  It  is  confined  to  the  cathode  gas  diffusion  and  catalyst  layers 
in  contact  with  the  membrane  and  with  a  gas  distributor,  here 
shown  by  two  half  channels  and  a  shoulder.  Depending  on  the 
selected  boundary  conditions  at  the  channel  interfaces  the  oper- 


membrane 


reaction 


y=H2 


y=H1 


diffusion  layer 


x=0 


x=L 


Fig.  1 .  Two-dimensional  model  domain. 


ation  with  conventional  straight  channel  and  with  interdigitated 
gas  distributors  can  be  modeled. 

It  has  been  demonstrated  that  the  main  mass  transport  lim¬ 
itations  in  PEM  fuel  cells  occur  within  the  electrode  (gas  dif¬ 
fusion  and  catalyst)  layers  [3].  For  the  case  of  a  conventional 
flow  field,  air  flows  parallel  to  the  gas  diffusion  layer.  The  pre¬ 
dominant  mechanism  for  oxygen  transport  within  the  electrode 
is  diffusion.  Liquid  water  moves  due  to  capillary  forces.  In  the 
case  of  interdigitated  flow  fields  the  gas  is  forced  to  flow  through 
the  diffusion  layer,  hence  producing  a  forced  convective  mass 
transfer.  Liquid  water  flows  due  to  both  capillary  forces  and  the 
shear  force  of  the  gas  stream. 

It  has  been  found  that  the  kinetics  of  the  oxygen  reduction 
reaction  (ORR)  at  the  cathode  is  at  least  five  orders  of  magnitude 
slower  than  the  hydrogen  oxidation  reaction  (HOR)  [24].  The 
cell  overpotential  on  the  anode  side  is  thus  neglected,  consider¬ 
ing  the  operation  with  pure  hydrogen.  The  catalyst  layer  at  the 
cathode  is  treated  as  a  thin  layer  where  oxygen  is  depleted  and 
water  and  heat  are  produced,  defining  source  and  sink  terms, 
respectively.  The  depletion  and  production  rates  depend  on  the 
local  current  density,  which  is  described  by  the  Tafel  equation 
[25].  The  local  current  density  is  again  a  function  of  the  local 
oxygen  concentration.  No  limitations  for  the  proton  transport  are 
considered  within  the  catalyst  layer.  The  membrane  is  assumed 
to  be  fully  hydrated  with  liquid  water  and  to  have  a  constant 
ionic  conductivity  for  the  determination  of  ohmic  losses  associ¬ 
ated  with  ionic  transport. 

For  the  description  of  the  transport  of  both  gas  and  liquid 
phases  in  the  porous  cathode  electrode  a  continuum  approach 
is  adopted,  see  e.g.,  [23].  It  is  also  assumed  that  the  system  is 
in  local  thermodynamic  equilibrium.  The  porous  electrode  is 
assumed  to  be  isotropic.  This  is  a  simplification  for  the  numer¬ 
ical  modeling.  The  authors  are  aware  that  the  material  can  be 
anisotropic  to  a  certain  degree  due  its  structure.  Furthermore,  an 
extended  Darcy’s  law  is  used  to  describe  the  movement  of  both 
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gas  and  liquid  phases  within  the  porous  system.  The  membrane 
is  considered  impermeable  for  the  gases.  The  ohmic  heating  in 
the  current  collector  plate  and  in  the  GDL  is  neglected  due  to 
their  high  conductivity.  A  part  of  the  heat  produced  in  the  elec¬ 
trochemical  reaction  is  transported  by  conduction  through  the 
membrane. 

2.2.  Governing  equations 

2.2.1.  Multiphase  flow  in  porous  media 

Using  an  Eulerian  approach,  we  formulate  a  balance  equation 
for  each  of  the  mass  components,  assuming  that  all  phases  are 
in  local  thermodynamic  equilibrium.  The  investigated  system  is 
composed  of  the  two  phases  (a)  liquid  water  and  gas,  represented 
by  ‘w’  and  ‘g’,  respectively,  and  the  three  components  (K)  H2O, 
O2  and  N2.  The  mass  balance  equation  per  component  can  be 
written  as 


d(J2a  (pSapaX^) 
3 1 


y.  di v{pax%va  +  JK]  -  qK  =  0, 


a  e  {w,  g},  K  e  {H20,  02,  N2}, 


(4) 


where  4>  represents  the  porosity,  Sa  the  phase  saturation,  defined 
as  the  volumetric  portion  of  the  pore  space  occupied  by  the  phase 
a,  pa  the  phase  density,  the  molar  fraction  of  component 
K  in  phase  a ,  va  the  Darcy  velocity,  JK  the  diffusive  flux  of 
component  K  and  qa  represents  a  source  (qa  >  0)  or  sink  (qa  < 
0)  term,  which  is  only  specified  within  the  reaction  layer. 

The  phase  velocities  are  determined  using  an  extended 
Darcy’s  law  for  multiphase  systems  [26]: 


va  = - K (grad pa  -  pa g), 

l^a 


(5) 


where  km  represents  the  relative  permeability  of  phase  a,  //„ 
the  dynamic  viscosity,  K  the  absolute  permeability  of  the  porous 
medium,  pa  the  pressure  of  phase  a  and  g  represents  the  gravity 
vector.  The  gravitation  effect  is  neglected  by  assuming  that  the 
gas  flow  in  the  channel  is  downwards,  thus  the  model  domain  of 
Fig.  1  corresponds  to  a  horizontal  plane. 

The  capillary  pressure  represents  the  discontinuity  in  pressure 
that  exists  across  the  interface  between  the  two  phases  when  they 
are  in  contact  in  the  interstices  of  a  porous  medium: 


Pc  =  Pg  ~  P w, 


(6) 


where  pg  is  the  pressure  in  the  gas  phase  and  pw  is  the  pressure  in 
the  water  phase.  Within  a  capillary  tube  of  radius  r,  the  Laplace 
equation  describes  the  capillary  pressure  as 


Pc  = 


2 a  cos  6 
r 


(7) 


where  a  is  the  surface  tension  and  6  is  the  contact  angle.  By 
convention,  9  (0°  <  9  <  180°)  is  measured  through  the  denser 
fluid.  When  9  <  90° ,  the  fluid  is  said  to  wet  the  solid  and  is  called 
a  wetting  fluid.  When  0  >  90°,  the  fluid  is  called  a  non-wetting 
fluid. 

In  the  studied  system,  with  a  hydrophobic  gas  diffusion  ma¬ 
terial,  the  liquid  water  phase  represents  the  non-wetting  phase 


with  9  >  90°  and  the  gas  is  the  wetting  phase,  hence  the  liq¬ 
uid  pressure  is  larger  than  the  gas-phase  pressure.  The  capillary 
pressure  is  thus  negative. 

Eq.  (7)  also  elucidates  that,  in  the  case  of  capillary  pres¬ 
sure  effects  in  hydrophobic  media,  the  wetting  fluid,  i.e.,  the  gas 
phase,  will  occupy  the  smaller  pores,  whereas  the  non-wetting 
fluid,  i.e.,  the  water  phase,  will  occupy  pores  with  a  larger  radius. 

In  a  macroscopic  description  of  porous  media  a  relationship 
pc  =  pc(Sw)  is  needed.  Hydrophilic  media  have  been  inten¬ 
sively  investigated  in  problems  of  earth  sciences  and  environ¬ 
mental  engineering  in  the  last  four  decades.  There  are  many 
models  which  have  been  obtained  from  experiments  with  dif¬ 
ferent  soil-like  (sandstone)  materials.  Most  commonly  used  are 
the  Leverett  approach  [27],  the  Van  Genuchten  model  [28]  and 
the  Brooks  and  Corey  model  [29].  They  include  fitting  param¬ 
eters  like  a  pore-grain-size-form  factor  and  Leverett  [27]  also 
includes  a  temperature  dependence  through  the  surface  tension 
and  contact  angle. 

Hydrophobic  media  have  been  considerably  less  investigated. 
[20,21]  used  the  Leverett  approach  and  included  a  contact  angle 
>90°  to  obtain  a  capillary  pressure-saturation  relationship  for  a 
hydrophobic  medium.  In  our  work,  we  follow  another  approach. 
Based  on  the  concept  of  capillary  pressure  at  the  microscale  (Eq. 
(7)),  which  is  widely  applied  to  determine  the  pore  size  distribu¬ 
tion  in  any  porous  medium  (e.g.,  mercury  intrusion  porosimetry) 
we  obtain  the  capillary  pressure-saturation  data  that  is  generated 
directly  during  the  measurement  and  correct  it  for  a  water-air 
system  considering  the  differences  in  surface  tension  and  contact 
angle.  In  this  way,  one  obtains  a  set  of  data  which  is  compara¬ 
ble  with  the  behaviour  of  water  in  the  porous  carbon  material. 
See  Section  4  for  a  detailed  description  and  the  results  of  this 
method. 

The  relative  permeability  of  a  phase  describes  how  the 
presence  of  one  phase  influences  or  disturbs  the  flow  behaviour 
of  the  other  phase  and  vice  versa.  In  the  earth  sciences  literature 
this  relationship  is  commonly  derived  from  the  capillary 
pressure-saturation  curve  which  fits  the  hydrophilic  material 
data.  The  most  common  form  of  the  relative  permeability  is 
represented  by  a  cubic  function  of  the  saturation.  The  function 
krw(Sw)  for  the  liquid  water  phase  (in  a  sandstone  the  wetting 
phase)  is  characterized  by  a  slight  increase  for  low  saturations 
and  a  strong  increase  for  higher  saturations.  This  is  due  to  the 
fact  that,  in  the  case  of  low  saturations,  the  liquid  (wetting 
fluid)  fills  only  very  small  pores,  where  flow  movements  are 
nearly  impossible  because  of  the  strong  molecular  attraction, 
whereas  the  largest  pores  are  only  filled  in  the  case  of  higher 
saturations. 

Because  in  hydrophobic  media  the  liquid  water  phase  fills 
first  the  larger  pores,  it  is  expected  to  obtain  a  faster  increase 
of  the  relative  permeability  of  water  at  low  saturations  than  by 
hydrophilic  media.  The  reason  is  that,  for  lower  saturations,  the 
larger  pores  are  filled  first,  whereas,  in  the  case  of  almost  full 
saturation,  the  remaining  very  small  pores,  which  are  filled,  have 
no  significant  influence  on  the  flow  behaviour.  For  this  reason 
it  was  decided  to  derive  the  relative  permeability-saturation  re¬ 
lationship  from  the  obtained  capillary  pressure-saturation  data 
for  the  present  system,  as  described  in  Section  4. 
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For  the  description  of  the  diffusive  component  transport  in 
each  phase  Fick’s  law  is  used.  For  the  gas  phase  it  can  be  used 
to  describe  the  diffusion  of  a  trace  component  in  a  carrier  gas, 
when  the  interdependencies  of  the  individual  mass  fluxes  of  the 
components  can  be  neglected.  Due  to  the  simplicity  of  Fick’s  law 
in  describing  explicitly  the  molar  diffusive  flux  of  each  compo¬ 
nent,  it  was  decided  to  make  this  assumption  and  to  insert  the 
flux  term  directly  into  the  mass  balances.  The  diffusive  flux  of 
each  component  is  thus 

JK  =  ~Dpmpa  gmdxK .  (8) 

Dpm  is  the  effective  macroscopic  diffusion  coefficient  for  com¬ 
ponent  K  in  the  porous  medium.  The  effective  diffusion  coeffi¬ 
cient  depends  on  the  tortuosity  r,  for  a  multiphase  system  it  also 
depends  on  the  gas  saturation: 

DU  =  \h°K  (9) 

where  DK  is  the  binary  diffusion  coefficient  for  the  component 
K  in  the  carrier  medium. 

A  single  energy  equation  is  formulated  for  the  fluid-filled 
porous  medium,  based  on  the  assumption  of  local  thermal  equi¬ 
librium  in  which  all  phases  have  the  same  temperature.  In  this 
way  heat  transport  terms  between  the  fluid  phases  can  be  ne¬ 
glected.  The  balance  equation  can  be  formulated  as  follows 

3(V„  Sapaua)  3pscsT 

V  3 1  ^  3 1 

-  div(kpmgradr)  -  ^  div{km—haK(gra.dpa  -  pa g)} 

a 

-  div f  DpmPghg  gradx^ }  -  qh  =  0,  (10) 

K 

where  ua  is  the  internal  energy  of  phase  a,  cs  is  the  heat  capacity 
of  the  solid  matrix,  kpm  is  an  average  thermal  conductivity  of 
the  fluid  filled  porous  medium,  ha  is  the  enthalpy  of  phase  a  and 
qh  represents  a  heat  source  or  sink. 

2.2.2.  Electrochemical  reaction 

The  cell  potential  can  be  determined  by 

t^cell  =  U0c  t]  a  Pc  IR spec;  (11) 

where  U0c  is  the  open  circuit  voltage  (here  the  reversible  cell 
voltage  is  used),  ;/a  and  t]c  the  anodic  and  cathodic  overvoltages, 
i  the  current  density  and  Rspec  is  the  specific  ohmic  resistance 
for  the  membrane,  electrodes  and  contacts  determined  experi¬ 
mentally  by  impedance  spectroscopy.  These  values  are  listed  in 
Table  3.  The  anodic  overvoltage  is  neglected  in  the  present  work, 
due  to  its  insignificant  contribution. 

The  transfer  overvoltage  at  the  cathode  trails  is  given  using 
the  Tafel  kinetics  for  the  oxygen  reduction  reaction  [25]: 


i  =  z'oef/v  exp 


-A Gt 
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(12) 


where  irf  represents  the  exchange  current  density  at  a  reference 
state,  /v  a  surface  increasing  factor  which  describes  the  rela¬ 
tion  between  the  catalytically  active  and  the  geometrical  surface 
area  of  the  electrode,  AGq  the  free  activation  enthalpy,  a  the 
load  transfer  factor,  z  the  number  of  transferred  electrons  per 
molecule  of  reactant  and  F  is  the  Faraday  constant.  The  factor 
(1  —  .Sw)  accounts  for  the  covering  effect  of  the  catalyst  sites  by 
the  liquid  phase  [15].  p° 2  is  the  oxygen  partial  pressure  in  the 
reaction  layer. 

Defining  q*  as 
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r  RT 

~ fv  +  ln(1  ~  Sw^ 


(14) 


Moreover,  for  the  space  discretization  it  will  be  considered  that 
the  open  circuit  voltage  at  the  reaction  place,  calculated  with  the 
local  concentration,  is  not  equal  to  the  open  circuit  voltage  at 
cell  inlet,  calculated  for  the  cell  inlet  conditions.  This  difference, 
considered  as  an  equilibrium  overvoltage  qu  [30],  is  given  by 


(15) 


where  t>;  are  the  stoichiometric  coefficients  of  the  reacting 
species  in  Eq.  2.  Now  the  actual  cathodic  overpotential  is 


tic  —  r/tians  +  Pequ- 


(16) 


Inserting  Eq.  14  and  15  in  Eq.  1 1  results: 
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For  a  multidimensional  case  the  local  current  density  i  for  a  given 
constant  cell  voltage  results  in  an  implicit  relation: 


R, 
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The  oxygen  consumption  rate  in  the  electrochemical  reaction  is 
given  by  Faraday’s  law: 


a£2  = 


i 

4 F’ 


(19) 


It  is  implemented  as  a  (volumetric)  sink  term  in  the  reaction 
layer.  In  the  same  way,  the  water  production  rate  is  given  by 


Nr  -  — 
^vh2o  -  -  -  • 


IF 


(20) 
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In  addition  the  water  flux  through  the  membrane  due  to  electro- 
osmotic  drag  is  being  considered.  A  water  transport  number 
(?h2o)  f°r  a  fully  hydrated  membrane  has  been  assumed  (see 
Table  3): 

<0  =  bh0^.  (21) 

The  heat  production  by  the  reaction,  considered  as  a  heat  source, 
is  calculated  as 

qR  =  (Uth  -  Uce\\)i.  (22) 

It  is  assumed  that  a  part  of  the  heat  generated  in  the  cathode 
reaction  layer  flows  through  the  membrane  to  the  anode  side, 
depending  on  the  thermal  conductivity  of  the  membrane. 

2.3.  Boundary  conditions 

Assuming  thermodynamic  equilibrium  and  according  to  the 
Gibbsian  phase  rule  this  non-isothermal  system  with  two  phases 
and  three  components  requires  four  independent  primary  vari¬ 
ables  to  be  uniquely  described.  Depending  on  the  phase  state 
of  the  system  at  a  local  point,  some  variables  may  need  to  be 
switched  (see  Section  3.2).  In  this  section  the  initial  and  bound¬ 
ary  conditions  are  explained  assuming  that  both  phases  (liquid 
and  gas)  are  present. 

The  boundary  conditions  available  in  the  numerical  simulator 
MUFTE  _UG  for  the  Box  and  CVFE  discretization  methods  (see 
Section  3)  are  of  types  Dirichlet  and  Neumann  (see,  e.g.,  [26]). 

The  governing  equations  are  identical  for  both  conventional 
and  interdigitated  flow  fields.  The  difference  exists  in  the  bound¬ 
ary  conditions  for  the  primary  variables  at  the  interface  between 
the  gas  channel  and  the  diffusion  layer.  In  all  other  boundaries, 
the  two  types  of  flow  fields  have  the  same  boundary  conditions. 

For  0  <  x  <  L  and  y  =  0  or  v  =  //,  symmetry  conditions 
are  used  because  y  —  0  and  H  correspond  to  the  center  of  the 
flow  channels  (see  Fig.  1).  Neumann  zero  boundary  conditions 
are  applied  (the  fluxes  of  all  mass  components  and  the  enthalpy 
flux  are  set  equal  to  zero). 

At  the  interface  with  the  membrane  (0  <  y  <  H  and  x  =  0) 
also  Neumann  zero  boundary  conditions  are  taken,  because  the 
electrochemical  reaction  is  assumed  to  occur  within  the  reaction 
layer  thickness.  Also  the  electro-osmotic  drag  flow  of  water  is 
for  simplicity  assumed  to  occur  in  the  reaction  layer,  which  is 
described  by  a  source  term. 

At  the  interface  in  contact  with  the  shoulder  of  the  current 
collector  (HI  <y<  H2  and  x  —  L)  also  Neumann  zero  con¬ 
ditions  are  applied  for  all  mass  components.  For  the  heat  term, 
a  constant  heat  conduction  flow  is  assumed  to  occur  through  the 
shoulder,  calculated  by 

■\  shoulder 

p^A7-iing),  (23) 

where  xshoulder  and  5shoulder  are  the  thermal  conductivity  and  the 
thickness  of  the  shoulder,  respectively. 


2.3.1.  Conventional  gas  distributor 

For  H2  <  y  <  H  and  0  <  y  <  HI  and  x  —  L  (i.e.,  at  the 
inlet  of  the  reacting  gases)  Dirichlet  boundary  conditions  for 
oxygen  concentration  in  the  gas  phase,  liquid  water  saturation, 
total  gas  pressure  and  temperature  are  set: 

x°2  =  x°2'inlet,  Sw  =  o,  pg  =  pfe\  T  =  rnlet  (24) 

In  this  work,  treating  hydrophobic  materials  it  is  for  simplicity 
assumed  that  the  liquid  water  saturation  at  the  channel  interface 
is  near  to  zero.  In  the  real  operation,  depending  on  the  actual 
wettability  of  the  gas  diffusion  material  and  on  the  water  vapor 
content  of  the  gas  stream,  a  specific  value  is  expected  to  be  es¬ 
tablished.  One  approach  to  estimate  this  saturation  value  at  the 
interface  could  be  obtained  from  a  local  equilibrium  assumption, 
based  on  sorption/desorption  isotherms  determined  experimen¬ 
tally  for  these  materials  under  different  operating  conditions  of 
water  vapor  saturation  in  the  gas  phase.  This  is  out  of  the  scope  of 
the  present  work  and  should  be  studied  in  further  investigations. 

2.3.2.  Interdigitated  gas  distributor 

For  H2  <  y  <  H  and  x  —  L  (flow  inlet)  Neumann  fluxes  of 
all  mass  components  are  set.  For  the  energy  equation  a  Dirichlet 
condition  for  the  inlet  temperature  is  used. 

For  0  <  y  <  HI  and  x  —  L  (flow  outlet)  the  fluid  phases 
can  break  through  into  the  gas  channel  at  the  outlet  pressure 
/7°utlet.  The  difficulty  in  the  modeling  arises  from  the  fact  that 
this  boundary  cannot  be  described  correctly  by  either  a  Dirichlet 
or  a  Neumann  boundary  condition.  Instead,  we  had  to  choose 
a  modified  model  domain  that  allows  a  quasi-simulation  of  the 
conditions  at  the  outlet  in  the  gas  channel.  We  extended  the 
domain  in  the  horizontal  direction  as  shown  in  Fig.  2.  At  the 
new  boundary  (0  <  y  <  HI  and  x  =  /.extension)  a  liquid  water 
saturation  of  zero  was  assumed.  We  assigned  to  this  extended 
portion  high  flow  properties  (permeability)  in  order  to  make  sure 
that  the  outlet  of  the  diffusion  layer  is  influenced  only  minimally 
by  the  righter  Dirichlet  condition  for  the  water  saturation.  For 
the  oxygen  transport,  it  is  assumed  that  the  gas  is  completely 
mixed  in  the  outlet  channel  at  the  interface  with  the  gas  diffusion 
layer  (0  <y<  H 1  and  x  =  L ).  At  the  outlet  of  the  extended 
portion  (0  <y<  HI  and  x  =  /.extension)  the  diffusion  term  is 
eliminated,  to  make  sure  that  back-diffusion  has  no  influence  on 
the  real  model  domain. 

For  the  energy  equation,  a  free  boundary  for  the  enthalpy 
flow  is  chosen  at  (0  <  y  <  HI  and  x  =  L),  so  that  the  heat  flux 
is  allowed  out  of  the  domain  as  it  occurs  due  to  the  respective 
gradient  given  at  the  boundary. 

3.  Numerical  solution  technique 

The  system  of  flow  and  transport  equations  as  explained  in 
the  previous  sections  puts  a  high  demand  on  the  numerical  so¬ 
lution  techniques.  The  resulting  system  of  partial  differential 
equations  is  highly  coupled  and  non-linear.  The  equations  are  of 
a  mixed  hyperbolic  and  parabolic  character  with  a  varying  dom¬ 
inance  depending  on  the  influences  of  capillary  pressure  and 
the  ratio  of  convective  to  diffusive  fluxes.  The  equations  have 
been  implemented  into  the  numerical  simulator  MUFTE_UG 
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Fig.  2.  Extended  two-dimensional  model  domain. 

[31]  which  was  developed  as  a  research  code  for  multiphase 
flow  in  porous  media.  The  model  domain  shown  in  Fig.  1  is, 
after  multigrid  refinement,  divided  into  6400  equally  sized  rect¬ 
angular  grid  cells.  There  are  40  divisions  in  .r-direction  and  160 
divisions  in  y-direction. 

3.1.  The  numerical  simulator  MUFTEJJG 

MUFTE_UG  is  a  sophisticated  numerical  tool  that  is 
continuously  developed  in  a  joint-venture  project  of  the 
Interdisciplinary  Center  for  Scientific  Computing,  University  of 
Fleidelberg,  (part:  UG)  and  the  Institute  of  Hydraulic  Research, 
University  of  Stuttgart  (part:  MUFTE).  The  toolbox  UG, 
standing  for  Unstructured  Grids,  is  the  developing  platform.  It 
provides  a  number  of  powerful  numerical  algorithms  for  solving 
systems  of  partial  differential  equations  on  unstructured  grids, 
cf.  [32]  or  [33].  MUFTE  stands  for  multiphase  flow  transport 
and  energy  model.  It  can  be  applied  to  simulate  isothermal 
and  non-isothermal  multiphase  multicomponent  flow  in  porous 
and  heterogeneous  or  fractured  media.  MUFTE  contains  the 
problem- specific  implementation  of  the  equations,  including 
discretization  methods,  constitutive  relationships,  and  bound¬ 
ary  conditions.  Both  UG  and  MUFTE  have  strictly  modular 
structure  allowing  the  combination  of  different  problem  classes 
with  various  numerical  solution  methods.  A  good  overview 
of  the  available  model  concepts  and  algorithms  is  given,  for 
example,  by  Helmig  [26]  or  Class  et  al.  [34]. 

3.2.  Adaptive  choice  of  the  primary  variables 

The  non-isothermal  system  of  the  two  phases  gas  and  water 
including  the  three  mass  components  water,  oxygen  and  nitro¬ 
gen  is  described  by  the  mass  balance  equations  and  the  thermal 
energy  equation  as  explained  in  Section  2.2.1.  Depending  on 
the  number  of  phases  locally  present  in  a  control  volume  differ¬ 
ent  phase  states  are  distinguished.  For  this  work,  two  relevant 
phase  states  are  the  state  ‘both  phases’  and  the  state  with  only 


the  gas  phase  present.  If  both  phases  exist,  water  saturation  Sw, 
gas-phase  pressure  p„,  temperature  T,  and  oxygen  mole  frac¬ 
tion  in  the  gas-phase  ;cg 2  are  chosen  as  the  primary  variables.  A 
single-phase  gas  state  requires  a  switching  of  the  primary  vari¬ 
able  saturation  to  an  additional  gas-phase  mole  fraction.  The 
switching  algorithm  and  further  details  can  be  found  in  [34]. 

3.3.  Box  discretization 

Eqs.  (4)  and  (10)  are  spatially  discretized  with  the  so-called 
Box  method,  which  represents  a  strictly  mass-conservative  Sub- 
domain  Collocation  Method,  cf.  [26].  The  Box  method  requires 
the  construction  of  a  secondary  mesh.  This  is  achieved  by  con¬ 
necting  the  centers  of  gravity  of  the  elements  with  the  mid-points 
of  the  element  edges.  By  this  each  node  is  assigned  an  unique 
control  volume.  On  the  other  hand,  each  element  contains  a 
number  of  sub-control  volumes  equal  to  the  number  of  nodes. 
The  discretization  scheme  can  be  derived  using  the  principle  of 
weighted  residuals  applied  on  the  primary  finite  element  mesh 
with  piece- wise  constant  weighting  functions  for  the  control  vol¬ 
umes  (boxes)  on  the  secondary  mesh.  We  apply  a  mass  lumping 
for  avoiding  non-physical  oscillations  of  the  solution.  The  ad- 
vective  terms  in  the  equations  are  computed  with  a  fully  upwind 
technique.  This  implies  an  upstream  evaluation  of  mobilities, 
mole  fractions  and  enthalpies. 

3.4.  Time-step  procedure  and  solution  methods 

The  system  of  partial  differential  equations  resulting  from 
Eqs.  (4)  and  (10)  is  highly  non-linear.  It  is  linearized  using  a 
damped  in-exact  Newton-Raphson  method,  as  described  in  [35]. 

For  the  time  discretization  we  use  a  fully  implicit  Euler 
scheme  applied  to  the  storage  terms.  We  use  here  a  stabilized 
bi-conjugated  gradient  method  (Bi-CGStab)  optionally  with  an 
extended  multigrid  preconditioner,  cf.  [33,34]  as  linear  solver 
for  the  Jacobian  system. 

4.  Parameter  determination 

Properties  of  the  commercial  GDL  material  Double  Sided 
ELAT  from  E-TEK,  Inc.,  needed  in  the  model  are  determined 
experimentally.  These  are  described  in  the  following.  For  illus¬ 
tration  purposes  some  SEM  photographs  of  the  top  and  a  cross 
section  of  an  ELAT  electrode  are  shown  in  Figs.  3-5,  respec¬ 
tively. 

4.1.  Contact  angle 

The  wetting  properties  of  the  gas  diffusion  material  ELAT- 
DS  were  investigated  by  measuring  the  contact  angle  with  the 
sessile  drop  method  [36,37].  Two  different  samples  of  the  ma¬ 
terial  from  the  same  production  batch  are  cut  out  with  a  size  of 
50  mm  x  50  mm.  The  contact  angle  is  measured  at  three  differ¬ 
ent  positions  on  each  sample.  To  measure  the  contact  angle  a 
droplet  of  deionised  water  with  3  mm  diameter  is  placed  on  the 
surface  of  the  sample.  Then  a  picture  from  the  side  is  taken  with 
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Fig.  3.  SEM  micrograph  of  an  ELAT-DS  electrode  -  top  view  -  50  x  magnifi-  Fig.  5.  SEM  micrograph  of  an  ELAT-DS  electrode  -  cross-section  -  lOOx 
cation.  magnification. 


a  CCD  camera  [36].  The  image  is  analyzed  using  the  Software 
DSA1  v.  1.80  Drop  Shape  Analysis  from  Kriiss  GmbH,  which 
calculates  the  contact  angle  by  fitting  a  tangent  on  the  drop  pro¬ 
file  at  the  three-phase  point  where  the  liquid  surface  contacts 
the  solid  surface.  For  the  ELAT-DS  gas  diffusion  layer  the  mea¬ 
sured  contact  angle  for  water  is  143°  with  a  standard  deviation 
of  2°.  For  mercury  the  contact  angle  has  a  value  of  141°  on  the 
ELAT-DS  electrode. 

4.2.  Porosity 

The  pore  size  distribution  of  the  investigated  gas  diffusion 
material  ELAT-DS  is  measured  by  mercury  intrusion  porosime- 
try  [38]  using  a  Pascal  240  Porosimeter  from  Porotec  GmbH. 
This  method  forces  mercury  with  its  high  surface  tension  to  pen¬ 
etrate  the  pores  of  the  gas  diffusion  layer.  Setting  the  amount  of 
mercury  uptake  as  a  function  of  pressure,  allows  for  the  determi- 


Fig.  4.  SEM  micrograph  of  an  ELAT-DS  electrode  -  top  view  -  10  000  x  mag¬ 
nification. 


nation  of  the  pore  size  distribution.  The  results  of  the  measured 
pore  size  distribution  are  shown  in  Fig.  6. 

To  calculate  the  porosity,  the  thickness  of  the  sample  is  de¬ 
termined  via  SEM.  For  ELAT-DS  a  thickness  of  0.45  mm  was 
obtained.  The  resulting  porosity  of  ELAT-DS  is  78%. 

4.3.  Tortuosity 

The  Wicke-Kallenbach  diffusion  cell  [39]  allows  for  the 
investigation  of  the  diffusivity  of  dry  gases  through  the  gas 
diffusion  material.  The  diffusion  cell  consists  of  two  chambers, 
which  are  separated  by  the  investigated  gas  diffusion  layer 
ELAT-DS.  The  pressure  difference  between  the  chambers  is 
set  to  zero.  Well-defined  gas  flows  through  each  chamber  pass 
the  surface  of  the  investigated  gas  diffusion  layer.  With  a  gas 
chromatograph  each  gas  flow  composition  is  measured  at  the 
inlet  and  outlet  of  the  diffusion  cell  to  determine  the  diffusive 
fluxes  through  the  gas  diffusion  layer.  The  effective  diffusivity 
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is  related  to  the  binary  diffusion  coefficient  by 

«eff=  -DtJ.  (25) 

r 

Combining  /J)e ff  from  the  Wicke-Kallenbach  experiment  and 
the  porosity  <j>  measured  by  mercury  intrusion  porosimetry, 
an  experimentally  determined  value  for  the  tortuosity  r  is 
obtained  from  the  above  equation.  The  tortuosity  comprises 
all  factors  inhibiting  species  transport  within  the  dry  porous 
diffusion  layer,  such  as  indirections  and  dead-end  pores.  The 
measured  effective  diffusion  coefficient  of  H2  in  nitrogen  is 
5.533  x  10“6±  1.019  x  10“6m2s-'.  The  effective  diffu¬ 
sion  coefficient  for  O2  in  nitrogen  is  1.91  x  10-6  ±  0.68  x 
10_6m2s_l.  The  errors  in  the  effective  diffusion  coefficient 
are  mainly  due  to  the  quality  of  the  calibration  gases  for  the 
gas  chromatograph  (uncertainty  in  the  exact  concentration). 
Therefore,  the  calculated  tortuosity  of  ELAT-DS  is  3  ±  1. 

4.4.  Capillary  pressure 

For  the  determination  of  the  capillary  pressure-saturation  re¬ 
lationship  the  mercury  injection  porosimetry  method  was  also 
applied.  Mercury  is  normally  a  non-wetting  fluid.  The  core  is 
placed  in  a  chamber,  which  is  then  evacuated  to  a  very  low  pres¬ 
sure,  and  mercury  is  forced  into  the  core  under  pressure,  as  al¬ 
ready  mentioned  in  Section  4.2.  The  volume  of  mercury  injected 
at  each  pressure  determines  the  non-wetting  fluid  saturation  (tak¬ 
ing  into  account  the  volume  of  the  cell  and  the  bulk  volume  of  the 
sample).  The  rest  gas  represents  the  wetting  phase.  The  proce¬ 
dure  is  continued  until  the  core  sample  is  saturated  with  mercury, 
or  pressure  reaches  a  predetermined  value.  Applying  a  correc¬ 
tion  permits  to  obtain  the  values  for  a  water-air  system,  taking 
into  account  the  different  surface  tension  and  contact  angle,  as 
proposed  in  [38,40]: 

Pc  water— air  °water  cos  $water 

-  =  - .  (20) 

Pc  Hg— air  CTHg  cos  0\  jg 

Fig.  7  shows  the  obtained  corrected  data  for  a  water-air  system 
for  the  ELAT-DS  sample,  with  a  zoom  for  the  region  of  low 
capillary  pressures. 

The  hysteresis  typical  for  the  drainage  and  imbibition  pro¬ 
cesses  is  clearly  shown  in  Fig.  7.  As  one  can  observe,  during 
the  drainage  process  the  liquid  is  held  back  in  the  pores  of  the 
sample,  yielding  a  higher  liquid  phase  saturation.  At  the  end  of 
the  drainage  measurement,  a  liquid  amount  is  still  present  in  the 
sample,  which  represents  the  residual  saturation. 

It  should  be  mentioned  that  the  terms  “imbibition”  and 
“drainage”  will  be  used  throughout  this  contribution  in  their 
colloquial  meaning,  namely  imbibition  for  water  uptake  and 
drainage  for  de-watering.  This  is  in  contrast  to  literature  where 
drainage  is  often  defined  as  the  replacement  of  the  wetting 
fluid  (here  air)  through  the  non-wetting  fluid  (here  water).  Since 
the  electrode  materials  studied  have  both  hydrophobic  and  hy¬ 
drophilic  properties  which  is  discussed  below,  the  literature  def¬ 
inition  would  lead  to  confusion  if  the  material  properties  change 
between  wetting  and  non-wetting.  Given  the  capillary  pressure 
vs.  saturation  curve  of  Fig.  7,  two  points  need  to  be  discussed: 


water  phase  saturation  Sw  [-] 

0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1 


Fig.  7.  Capillary  pressure-saturation  data  for  DS-ELAT  electrode  corrected  for 
a  water-air  system.  The  insert  shows  the  relevant  range  of  capillary  pressure  and 
a  fitting  for  the  imbibition  curve. 


( 1 )  Which  curve,  imbibition  or  drainage  from  Fig.  7  should  be 
used  for  the  simulation? 

(2)  It  is  well  known  that  any  electrode  material  is  not  com¬ 
pletely  hydrophobic  but  takes  up  a  certain  amount  of  water 
at  zero  (capillary)  pressure,  whereas  mercury  is  completely 
unwetting.  The  conversion  of  the  mercury  porosimetry  data 
to  water  using  Eq.  (7)  does  not  take  this  difference  into 
account.  The  imbibition  or  drainage  curves  for  water  will 
therefore  in  reality  not  level  off  at  zero  saturation  as  in  Fig.  7 
but  extend  into  positive  capillary  pressures  as  in  Fig.  8. 

Since  at  present  no  detailed  experimental  studies  of  these 
effects  are  available,  we  will  in  the  following  consider  two  lim¬ 
iting  cases.  The  first  case  is  the  completely  unwetting  behaviour, 
characterized  by  the  imbibition  curve  fitted  to  Fig.  7.  The  sec¬ 
ond  case  is  the  (more  realistic)  case  of  partial  wetting,  assuming 
the  fit  to  the  drainage  curve  of  Fig.  8  with  an  extension  at  low 
saturation  into  positive  capillary  pressure. 

For  the  parametrization  of  the  imbibition  curve  the  following 
equation  has  been  used 

pc  =  Ai  exp(5iSw  +  Cj)  +  £>i(l  -  Sw).  (27) 

For  the  drainage  process,  further  measurements  with  water  have 
to  be  made  in  order  to  obtain  the  missing  data  at  low  saturations. 
In  this  work,  following  equation  has  been  used: 

pc  =  A2  exp(B2Sw  +  C2)  +  D2(l  -  Sw)  +  ^ ,  (28) 

Ow 

which  is  shown  in  Fig.  8  together  with  the  experimental  values. 

4.5.  Permeability 

4.5.1.  Absolute  permeability 

The  absolute  or  intrinsic  permeability  is  a  parameter  char¬ 
acterizing  only  the  porous  medium,  it  is  independent  of  the 
fluid  or  fluids  present.  For  the  determination  of  the  absolute 
through-plane  permeability  of  the  porous  electrode  material  an 
annulus-shaped  GDL  sample  was  situated  between  two  plates 
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Fig.  8.  Capillary  pressure-saturation  drainage  data  and  fitting  curve  for  DS- 
ELAT  electrode. 


of  a  metallic  cell.  Dry  air  was  fed  to  the  cell  at  a  controlled  rate, 
and  the  pressure  loss  between  gas  inlet  and  outlet  was  measured 
with  a  manometer  (Bexbill  England  FCO  12-Micromanometer). 
The  through-plane  permeability  for  the  investigated  ELAT-DS 
material  was  determined  as  5.20  x  10“ 11  m2,  which  is  a  lit¬ 
tle  lower  than  the  permeability  for  a  sandstone  (approximately 
x  10“ 10  m2)  but  is  in  the  range  of  other  commercial  GDL  mate¬ 
rials  as  reported  by  Ihonen  et  al.  [22] . 


4.5.2.  Relative  permeability 

As  mentioned  in  Section  2.2.1  the  relative  permeability  ac¬ 
counts  for  the  interaction  of  one  phase  on  the  flow  behaviour  of 
the  other.  In  a  continuum  approach  it  is  a  dimensionless  num¬ 
ber  between  0  and  1,  which  depends  on  the  phase  saturation. 
There  are  several  proposed  methods  for  the  determination  of  the 
relative  permeability-saturation  relationship  kIV/(Sw )  from  the 
capillary  pressure  data,  see  e.g.,  [41]  for  a  detailed  explanation 
of  the  different  models.  In  the  present  work  we  assumed  a  capil¬ 
lary  pore  model  of  parallel  tubes  with  constant  cross-section  in 
flow  direction  and  the  Burdine  [42]  approach  for  the  tortuosity 
dependency  on  the  saturation,  resulting  in 


krw  — 


kIg  — 


Sw  -  Sw 


1  - 


r\  f  Sw  d^w 

JO  „2 


PI 


f  1  d.Sw 

Jo  p2 


f  1  4  3  w 
Pc 

f1  dSw 

Jo  pi 


(29) 


(30) 


SWr  is  called  the  residual  saturation  of  the  water  phase.  It  repre¬ 
sents  the  amount  of  fluid  under  which  no  phase  continuity  exists 
anymore  and  the  phase  is  considered  immovable.  Fig.  9  shows 
the  obtained  curves  for  the  relative  permeability  of  each  phase 
using  the  imbibition  data. 

Note  that  krw  +  klg  <  1  due  to  the  pore  space  connectivity,  by 
which  the  flow  paths  are  interconnected  and  the  phases  interfere 
each  other. 


5.  Results  and  discussion 

The  model  was  used  to  evaluate  the  influence  of  the  flow 
configuration  and  capillary  properties  of  the  diffusion  layer 
material  on  the  cell  behaviour.  In  the  following  the  results  for 
the  imbibition  curve  are  shown  both  using  conventional  and 
interdigitated  type  gas  distributors.  Furthermore,  calculations 
using  the  extended  drainage  curve  from  Fig.  8  are  performed. 
It  should  be  noted  that  this  two-dimensional  model  is  limited 
to  the  channel  entrance  of  the  fuel  cell  and  can  not  take  into 
account  the  three-dimensional  effects  that  originate  along  the 
channels. 

5.1.  Imbibition  curve  results 

Tables  1-4  list  the  parameter  values  used  for  the  calculations. 
As  is  shown,  the  air  is  assumed  to  be  completely  humidified  with 
water  vapor  but  no  liquid  water  is  present  at  the  channel  interface. 
The  dimensions  of  the  channel  and  shoulder  heights  for  the  gas 
distributor  correspond  to  the  dimensions  used  in  our  fuel  cell 
with  an  active  area  of  25  cm2.  For  the  cell  voltage,  a  value  of 
0.5  V  is  used.  Other  conditions  and  parameters  used  in  the  study 
are  shown  in  these  tables  along  with  the  references  for  various 
parameters  used  in  the  model. 


Table  1 

Geometrical  and  gas  distributor  parameters 

Geometrical  parameters 

Length,  x  (L)  (cm) 

0.05 

Length,  v  (H)  (cm) 

0.2 

Gas  distributor  parameters 

Channel  width  (half)  (cm) 

0.05 

Shoulder  width  (cm) 

0.1 

Thermal  conductivity.  kshoulder  (WKT1  nr1) 

14.7 

Flow  field  plate  thickness,  ,5sholllder  (Cm) 

0.3 

A  T  for  cooling  (K) 

2 
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Table  2 

Electrode  parameters 

Diffusion  layer  parameters 

Diffusion  layer  thickness  (cm) 

0.045 

Porosity,  0 

0.78 

Absolute  permeability,  K  (  x  10~u) 

5.20 

Residual  water  saturation  for  imbibition 
process,  SWI 

0.05 

Residual  water  saturation  for  drainage 
process,  SWI 

0.41 

Contact  angle  water,  #water  (°) 

143 

Contact  angle  mercury,  #Hg  (°) 

141 

Surface  tension  of  water,  <rwate r  (N  m-1) 

0.0625 

Surface  tension  of  mercury,  <rHg  (N  m-1) 

0.480 

Constants  used  in  capillary  pressure-saturation 

Ai,BuC1,Dl  (Eq.  (27)) 

-1168.75,  8.5,  -0.2,  -700 

A2,  B2,  C2,  D2,  E2  (Eq.  (28)) 

-600,  25,-16,  -3300,  800 

Tortuosity,  r 

3 

Thermal  conductivity,  X.  [25] 

(wr'm'1) 

15.6 

Heat  capacity,  cp  [25]  (Jkg~*  KT1) 

710 

Reaction  layer  parameters 

Reaction  layer  thickness  [25]  (cm) 

0.005 

Porosity,  0  [25] 

0.07 

Tortuosity,  r  [25] 

5 

Table  3 

Parameters  for  the  electrochemical  reaction 

Reversible  voltage  at  T  =  70  °C,  I/rev  (V)  1.191 

Thennoneutral  voltage  at  T  =  70  °C,  17th  (V)  1.4836 

Free  activation  enthalpy  [43],  AGq  (  x  103  J mol-1)  73.0 

Cell  voltage  (V)  0.5 

Exchange  current  density  ORR  i^1  [43]  (x  10-8  A  cm'2)  1.87 

Reference  temperature  [43]  (K)  353.15 

Reference  partial  pressure  of  O2  [43]  (x  105  Pa)  5.0 

Transfer  coefficient  for  oxygen  reduction,  a  [25]  0.5 

Surface  increasing  factor,  /v  [25]  60 

Electro-osmotic  drag  coefficient,  tH->0  [44]  0.2327 

MEA  specific  resistance  Rs pec  [25]  (£2  nr)  0.25 

Membrane  thermal  conductivity,  XM  [45]  (WKT1  m_1)  0.43 

Membrane  thickness,  (p.m)  87.5 


5.1.1.  Conventional  flow  field 

For  the  conventional  flow  field,  the  gas  stream  flows  parallel 
to  the  diffusion  layer  and  the  oxygen  transport  mechanism  to 
the  reaction  layer  is  represented  by  diffusion.  Fig.  10  provides 
the  steady-state  concentration  profiles  of  oxygen,  liquid  water 
saturation  and  temperature  distribution  for  the  imbibition  curve. 
Fig.  1 1  shows  the  local  current  density  profile  along  the  electrode 
height  at  the  interface  membrane/catalyst  layer. 

At  the  channel/diffusion  layer  interface,  the  oxygen  concen¬ 
tration  is  set  to  inlet  boundary  condition  that  is  a  mole  fraction 
of  0.1775  for  oxygen  in  air  saturated  with  water  vapor  at  the 


Table  4 

Operational  conditions 

Inlet  gas  pressure  (x  105  Pa  (abs))  2.013 

Inlet  gas  temperature  (°C)  70 

Inlet  oxygen  mole  fraction  0. 1 775 

Inlet  water  vapor  mole  fraction  0. 1 546 

Inlet  liquid  water  saturation  0.01 

Inlet  stream  relative  humidity  (%)  100 


inlet  temperature  of  70  °C  and  pressure  of  2  bar.  From  this  inter¬ 
face  oxygen  diffuses  to  the  reactive  surface  through  the  porous 
gas  diffusion  layer  where  it  is  consumed  by  the  electrochemical 
reaction  to  produce  water.  Also  the  water  flow  coupled  to  the 
proton  transport  through  the  membrane  coming  from  the  anode 
side  is  accounted  for  at  the  reaction  layer.  The  water  accumu¬ 
lated  at  the  cathode  is  removed  in  general  by  two  mechanisms, 
i.e.,  capillary  flow  and  evaporation.  In  this  case,  the  air  flow  is 
completely  saturated  with  water  vapor,  so  that  capillary  action 
represents  the  dominant  mechanism. 

In  the  reactive  regions  under  the  channels,  the  reactants  and 
products  need  to  cover  only  the  thickness  of  the  diffusion  layer, 
therefore  mass  transport  is  fast.  Otherwise,  in  the  reactive  region 
under  the  shoulder,  the  reactants  and  products  have  to  diffuse 
across  the  thickness  and  height  of  the  diffusion  layer.  Hence, 
a  lower  availability  of  oxygen  and  a  greater  accumulation  of 
water  in  the  diffusion  layer  under  the  shoulder  is  expected  as 
can  clearly  be  seen  in  the  oxygen  mole  fraction  and  liquid  wa¬ 
ter  saturation  distribution  profiles  in  Fig.  10.  Since  the  current 
density  along  the  electrode  length  is  directly  affected  by  the 
oxygen  availability,  the  current  density  profile  shown  in  Fig.  1 1 
matches  the  oxygen  profile  in  the  catalyst  layer.  The  air  cathode 
average  current  density  for  the  whole  reaction  layer  for  the  im¬ 
bibition  curve  using  conventional  gas  distributors  is  calculated 
to  be  0.69  A  cm-2. 

The  liquid  water  saturation  in  the  cathodic  GDL  is  higher  in 
the  region  under  the  shoulder  and  increases  only  up  to  1 1.5  %. 
This  is  because  the  liquid  water  generated  under  the  shoulder 
accumulates  there  despite  the  lower  rate  of  water  generation 
in  this  region.  Due  to  the  longer  distance  to  the  channel  in  the 
region  under  the  shoulder,  the  liquid  water  can  not  be  removed 
as  quickly  as  in  the  region  directly  under  the  channels.  The 
liquid  water  flux,  induced  by  the  gradient  in  the  liquid  water 
saturation,  is  directed  from  the  catalyst  layer  into  the  channel, 
i.e.,  in  the  opposite  direction  to  the  gas  phase  (diffusion) 
flux. 

From  the  temperature  profile  (Fig.  10(c))  one  can  see  that  the 
temperature  is  also  higher  in  the  region  of  high  current  densities, 
that  is  the  reaction  zone  under  the  channels.  This  is  because  the 
shoulder  is  cooled  by  heat  conduction  into  the  current  collector. 
The  small  temperature  increase  within  the  GDL  depends  on  the 
thermal  conductivity  of  the  material  and  is  in  agreement  with 
other  simulation  studies,  e.g.,  [12]. 

5.1.2.  Interdigitated  flow  field 

When  an  interdigitated  gas  distributor  is  used,  oxygen  and 
water  are  transported  to  and  from  the  inner  layer  of  the  elec¬ 
trode  by  both,  diffusion  and  convection  induced  by  the  pressure 
drop  created  between  the  inlet  and  outlet  channels.  This  convec¬ 
tive  flow  through  the  electrode  reduces  the  gas  diffusion  distance 
to  and  from  the  catalyst  layer,  and  by  having  gas  flowing  under 
the  shoulders  of  the  gas  distributor,  the  electrode  active  area  un¬ 
der  the  shoulder  is  better  utilized  [13].  Also,  the  liquid  water 
removal  is  enhanced,  thereby  reducing  the  electrode  flooding 
problems.  These  benefits  help  to  extend  the  cell  operable  condi¬ 
tions  to  higher  current  densities  and  consequently  higher  power 
densities. 
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Fig.  10.  Oxygen  concentration,  Xg 2  (a),  liquid  water  saturation.  Sw  (b)  and  temperature.  T  (K)  (c)  distribution  in  the  electrode  for  the  imbibition  curve- 
gas  distributor  (X,  electrode  thickness  and  Y,  electrode  height). 
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Fig.  1 1 .  Current  density  distribution  at  the  catalyst  layer-membrane  interface 
for  the  imbibition  curve — conventional  gas  distributor. 


The  results  are  given  in  Figs.  12-14.  Like  for  the  conventional 
gas  distributor,  the  inlet  gas  stream  is  assumed  to  be  completely 
saturated  with  water  vapor,  but  no  liquid  water  is  present  in  the 
channel.  The  reactant  fluxes  were  taken  for  an  air  stoichiometric 
flow  rate  of  2  at  1  A  cm-2.  The  gas  inlet  temperature  is  set  to 
70  °C  and  the  gas  backpressure  is  set  to  2  bar. 


Fig.  12  shows  the  pressure  distribution  and  the  gas  veloc¬ 
ity  within  the  electrode.  At  the  inlet,  the  x-direction  velocity 
is  high  and  quickly  decreases  and  transfers  its  momentum  to 
the  y-direction  velocity.  The  same  occurs  at  the  outlet  where 
the  y-direction  velocity  decreases  and  transfers  the  momentum 
to  the  x-direction  velocity.  These  velocities  peak  near  the  in¬ 
let/shoulder  and  outlet/shoulder  corners  because  they  constitute 
the  shortest  distance  between  the  inlet  and  outlet  channels. 

Fig.  13  shows  the  oxygen  concentration,  liquid  water  satura¬ 
tion  and  temperature  distribution  profiles  within  the  electrode. 
Fig.  14  shows  the  local  current  density  profile  along  the  electrode 
height  at  the  interface  membrane/catalyst  layer. 

The  oxygen  mole  fraction  is  highest  in  the  inlet  region 
(0.0015  m  <  v  <  0.002  m).  From  there,  the  oxygen  concentra¬ 
tion  decreases  as  oxygen  moves  from  the  inlet  toward  the  cat¬ 
alyst  layer,  where  it  is  consumed  as  shown  in  Fig.  13(a).  The 
oxygen  mole  fraction  within  the  electrode  layer  continues  to 
decrease  along  the  y-direction  as  oxygen  is  continuously  con¬ 
sumed  along  the  catalyst  layer.  Water  is  generated  within  the 
catalyst  layer,  the  electro-osmotic  water  flow,  which  depends 
on  the  reaction  rate,  is  also  assumed  to  appear  within  the  thin 
catalyst  layer.  The  liquid  water  saturation  shown  in  Fig.  13(b) 
is  some  higher  in  the  region  under  the  shoulder  but  is  almost 
uniform  within  the  electrode.  The  maximum  liquid  water  sat¬ 
uration  for  this  case  is  about  8%  and  thus  a  little  lower  as  the 
value  obtained  for  the  conventional  gas  distributor.  The  current 
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Fig.  14.  Current  density  distribution  at  the  catalyst  layer-membrane  interface 
for  the  imbibition  curve — interdigitated  gas  distributor. 


density  profile.  Fig.  14,  corresponds  to  the  oxygen  profile  along 
the  interface  to  the  catalyst  layer  as  shown  in  Fig.  13(a).  The 
average  current  density  for  the  imbibition  curve  using  interdig¬ 
itated  gas  distributors  for  a  cell  voltage  of  0.5  mV  is  calculated 
to  be  0.75  A  cm-2. 

The  temperature  profile  (Fig.  13(c))  is  in  the  same  way  dif¬ 
ferent  from  the  conventional  gas  distributor.  In  this  case,  the 
gas  temperature  increases  throughout  the  electrode  due  to  the 
reaction  enthalpy  and  reaches  the  maximum  value  at  the  outlet 
in  spite  of  the  cooling  through  the  shoulder.  The  cooling  flow 
(considered  for  simplicity  as  constant)  is  for  this  case  not  high 
enough  to  reduce  the  temperature  of  the  outflowing  gas  stream. 
For  this  reason,  a  higher  temperature  increase  is  observed  for 
this  gas  distributor. 

When  comparing  the  results  for  both  gas  distributors  using 
the  imbibition  curve,  one  can  observe  a  lower  liquid  water 
saturation  within  the  electrode  when  the  interdigitated  flow 
field  is  used.  Due  to  the  faster  oxygen  transport  by  convection, 
the  oxygen  mole  fraction  and  thus  the  current  density  are  also 
higher  for  this  type  of  flow  distributor.  At  an  overvoltage  of 


Fig.  15.  Imbibition  curve  results  and  measured  polarization  curves  for  the  con¬ 
ventional  gas  distributor  operating  with  air. 


Fig.  16.  Imbibition  curve  results  and  measured  polarization  curves  for  the  in¬ 
terdigitated  gas  distributor  operating  with  air. 

200  mV  the  current  density  reaches  a  value  of  1.36  A  cm-2 
for  the  conventional  gas  distributor  and  1.82  A  cm-2  for  the 
interdigitated  one.  However,  e.g.,  Nguyen  et  al.  [13]  and  Kazim 
et  al.  [7]  report  much  higher  current  density  values  for  the 
interdigitated  gas  distributors.  This  is  due  to  the  fact  that  they 
assume  much  higher  values  for  the  pressure  drop  between  inlet 
and  outlet  channels,  as  Nguyen  et  al.  [13]  observed  from  their 
experiments. 

In  our  model,  at  the  inlet  boundary  we  assume  the  Neumann 
fluxes  for  the  reactant  species  so  that  the  pressure  drop  is  cal¬ 
culated.  The  obtained  value  of  about  160  Pa  is  much  lower  than 
the  value  reported  by  He  et  al.  [13]  and  Kazim  et  al.  [7],  for 
example,  that  is  in  the  range  of  700  Pa.  In  our  cell  laboratory 
measurements,  we  also  observe  higher  pressure  drops  than  pre¬ 
dicted  by  our  model,  but  this  pressure  drop  varies  significantly 
depending  on  the  humidity  of  the  gas  stream  and  the  current 
density,  i.e.,  the  amount  of  liquid  water  present  in  the  cathode, 
as  also  reported  by  He  et  al.  [46]  and  Thoben  and  Siebke  [47]. 
Because  the  liquid  water  saturation  distribution  is  a  variable  be¬ 
ing  calculated  for  the  GDL,  we  preferred  to  set  the  values  for  the 
fluxes  and  the  gas  backpressure  rather  than  the  pressure  drop  at 
the  electrode. 


Fig.  17.  Simulated  polarization  curve  for  a  partially  hydrophilic  electrode, 
drainage  curve,  and  ELAT-DS  measured  curve  for  comparison — conventional 
gas  distributor  operating  with  air. 
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The  low  value  obtained  for  the  pressure  drop  indicates  that 
in  practice  higher  liquid  water  saturations  are  to  be  expected  in 
the  GDL.  This  effect  can  not  be  reproduced  by  the  imbibition 
curve  in  Fig.  7.  Also  the  assumption  of  zero  water  saturation 
at  the  channel  interface  through  the  boundary  condition  has  to 
be  more  adequately  formulated  based  upon  measurements  of 
the  equilibrium  uptake  of  water  (water  saturation)  depending  on 
water  vapour  partial  pressure.  Moreover,  this  two-dimensional 
model  cannot  account  for  the  water  accumulation  along  the  chan¬ 
nels. 

5.2.  Comparison  with  experimental  data 

The  model  results  are  compared  to  a  set  of  experimental 
data  collected  on  a  25  cm2  cell  both  with  a  conventional  gas 
distributor  and  with  an  interdigitated  gas  distributor.  A  Nafion 
1135  membrane  (DuPont)  was  placed  between  two  electrodes 
of  type  Double  Sided  ELAT  from  E-TEK,  Inc.,  with  a  Pt 
loading  of  1 .0  mg  Pt  cm  2  (20%  Pt  on  Vulcan  XC-72)  and 
impregnated  with  Nafion  solution  with  a  loading  of  approx. 
1.5  mg  Nafion  cm-2.  Before  the  installation  in  the  test  facility, 
the  MEA  was  hot  pressed  [48],  The  polarization  curve  for 
this  MEA  was  generated  using  a  single-cell  assembly,  with  a 
conventional  flow  field  at  the  anode  side  and,  respectively,  an 
interdigitated  and  a  conventional  flow  field  for  the  cathode.  The 


cell  was  operated  at  70  °C,  the  anode  and  cathode  stream  were 
humidified  with  a  sparging  bottle  held  at  80  °C.  Hydrogen  was 
supplied  to  the  anode  and  air  to  the  cathode  at  a  stoichiometric 
flow  rate  of  1.5  and  2  at  1  A  cm-2,  respectively.  Both  gases  had 
a  pressure  of  2  bar. 

Figs.  15  and  16  show,  respectively,  the  experimental  results 
for  the  conventional  and  the  interdigitated  gas  distributor,  and 
the  model  results  using  the  imbibition  curve  from  Fig.  7  (In¬ 
sert).  In  both  diagrams  one  can  identify  two  different  effects.  At 
high  cell  potentials,  over  700  mV,  the  predicted  current  densities 
are  slightly  lower  than  the  measured  values.  The  reason  for  the 
deviations  may  be  inaccuracies  in  rate  constants  of  the  electro¬ 
chemical  reactions  in  this  kinetic  limited  region.  At  cell  voltages 
below  600  mV,  where  transport  limitations  become  significant, 
the  model  predictions  for  the  current  density  become  signifi¬ 
cantly  higher  than  the  experimental  values.  That  is  attributed  to 
the  low  liquid  water  saturation  values  predicted  by  the  model 
by  assuming  a  completely  hydrophobic  material.  In  practice, 
the  gas  diffusion  layer  material  is  only  partially  hydrophobic, 
and  can  also  turn  more  hydrophilic  under  cell  operation,  as  also 
reported  by  Reissner  et  al.  [49]. 

However,  as  shown  in  both  Figs.  15  and  16,  the  simulation 
results  predict  the  improvement  in  cell  performance  with  the 
interdigitated  gas  distributor  correctly,  showing  a  higher  limiting 
current  density  for  this  flow  configuration. 


(a)  X 


0.16 

0.14 


Fig.  18.  Liquid  water  saturation  and  oxygen  concentration  profile  in  the  electrode  at  0.5  V  cell  voltage  for  the  drainage  curve — conventional  gas  distributor  (X, 
electrode  thickness  and  Y,  electrode  height):  (a)  Sw;  (b)  Xg2. 
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Fig.  19.  Current  density  distribution  at  the  catalyst  layer-membrane  interface 
at  0.5  V  cell  voltage  for  the  drainage  curve — conventional  gas  distributor. 

5.3.  Drainage  curve  results 

In  order  to  obtain  the  system  response  by  allowing  higher 
liquid  saturations,  the  extension  of  the  drainage  curve  presented 
in  Fig.  8  is  used  for  the  calculations.  For  the  fitting  curve  selected, 


the  transition  from  the  hydrophobic  to  the  hydrophilic  region 
occurs  at  a  liquid  water  saturation  of  about  0.41.  This  value  is 
interpreted  as  the  lowest  water  saturation,  under  which  a  pressure 
in  the  gas  phase  has  to  be  applied  in  order  to  replace  the  water 
phase  trapped  in  the  pore  system. 

Fig.  17  shows  the  resulting  polarization  curve  obtained  for 
the  drainage  data,  compared  again  with  the  measured  data  for 
the  conventional  gas  distributor.  In  this  case  the  model  predic¬ 
tion  matches  the  experimental  values  very  good,  which  means 
that  in  practice  higher  liquid  water  saturations  are  to  be  ex¬ 
pected  in  the  electrode.  The  model  can  hence  correctly  re¬ 
produce  the  hindering  effect  of  the  liquid  phase  on  the  cell 
performance. 

The  drainage  curve  can  thus  better  represent  the  real  cell 
operation  than  the  imbibition  one.  This  is  also  expected  by  the 
fact  that,  after  a  given  operation  time,  liquid  water  will  be  present 
at  the  fuel  cell  cathode  and  the  drainage  process,  that  is  taking 
place  rather  than  purely  imbibition,  is  more  significant  for  the 
liquid  water  distribution. 

Fig.  18  shows  the  steady-state  oxygen  concentration  and 
liquid  water  saturation  profiles  considering  the  drainage  curve 
given  in  Fig.  8  for  a  cell  voltage  of  0.5  V.  As  expected,  high  liq¬ 
uid  water  saturation  values  are  obtained,  with  a  correspondingly 
lower  current  density  distribution  as  shown  in  Fig.  19,  with  an 
average  value  of  0.40  A  cm-2  at  0.5  V  cell  voltage. 
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Fig.  20.  Liquid  water  saturation  and  oxygen  concentration  profile  in  the  electrode  at  0.3  V  cell  voltage  for  the  drainage  curve — conventional  gas  distributor  ( X , 
electrode  thickness  and  Y,  electrode  height):  (a)  Sw;  (b)  Xg2. 
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Fig.  20  shows  the  steady-state  oxygen  concentration  and  liq¬ 
uid  water  saturation  profiles  considering  the  drainage  curve  for 
a  cell  voltage  of  0.3  V.  At  this  lower  voltage,  the  liquid  water 
saturation  in  the  electrode  increases  upon  60%. 

6.  Summary  and  conclusions 

A  two-dimensional,  two-phase,  non-isothermal,  multicom¬ 
ponent  model  was  developed  for  the  cathode  of  a  PEM  fuel 
cell,  which  can  be  applied  using  both  conventional  and  in- 
terdigitated  gas  distributors.  In  a  PEM  fuel  cell,  where  the 
membrane  needs  to  be  humidified  with  liquid  water  to  provide 
sufficient  protonic  conductivity,  the  water  management  at  the 
cathode,  where  water  is  being  produced  by  the  electrochem¬ 
ical  reaction  and  an  excess  can  hinder  the  oxygen  diffusion, 
becomes  a  difficult  task.  The  performance  of  the  fuel  cell  is 
highly  influenced  by  the  presence  of  liquid  water  and  its  removal 
rate,  particularly  at  high  current  densities.  The  water  accumu¬ 
lation  is  strongly  influenced  by  both  the  gas  flow  configuration 
and  the  wetting  properties  of  the  cathode  gas  diffusion  mate¬ 
rial. 

Due  to  the  lack  in  physical  properties  for  the  commonly  used 
electrodes,  the  most  important  properties  affecting  the  hydrody¬ 
namics  of  liquid  and  gas  flow  in  porous  media  have  been  de¬ 
termined  experimentally.  For  the  capillary  behaviour  an  indirect 
method  has  been  applied.  In  a  first  step  the  relative  permeabil¬ 
ity  data  have  been  derived  from  mercury  intrusion  porosimetry 
measurements  and  corrected  to  water  through  the  Laplace  equa¬ 
tion.  It  turned  out  however  that  a  more  detailed  consideration  of 
the  partial  wettability  of  the  electrode  backing  is  necessary  to 
obtain  agreement  between  simulation  and  experiments  at  high 
fuel  cell  loads.  In  future  work  the  partially  hydrophilic  behaviour 
of  the  gas  diffusion  material  should  be  closer  investigated  for  an 
air-water  system.  For  the  relative  permeability,  new  methods  for 
the  experimental  determination  are  going  to  be  part  of  the  work 
on  the  next  periods. 

The  increase  in  performance  due  to  the  faster  water  removal 
by  the  interdigitated  flow  distributor  is  correctly  predicted  by  the 
model.  However,  due  to  the  assumption  of  a  purely  hydrophobic 
character  of  the  gas  diffusion  material,  only  low  liquid  water 
saturations  are  obtained.  By  considering  a  partially  hydrophilic 
behaviour,  as  shown  during  the  drainage  process,  high  liquid  sat¬ 
urations  are  present  in  the  electrode  yielding  to  mass  transport 
limitations.  Thus,  the  model  predicts  a  decrease  in  performance 
due  to  the  stronger  hindering  effect  of  the  liquid  phase  on  the 
oxygen  transport  which  is  in  agreement  with  experimental  re¬ 
sults. 

The  suggested  model  predicts  correctly  the  tendency  of  the 
different  transport  processes  in  the  cathode,  whereas  the  main 
difficulties  in  describing  this  system  arise  from  the  determination 
of  the  most  important  flow  parameters,  in  particular  wettability, 
capillarity  and  relative  permeability,  which  will  be  closer  inves¬ 
tigated  for  these  materials  in  future  work.  On  the  other  side, 
also  the  correct  description  of  the  boundaries  to  the  gas  chan¬ 
nel  affecting  the  liquid  water  saturation  is  of  great  importance 
for  the  system  description.  The  liquid  water  saturation  value 
could  be  obtained  from  an  equilibrium  assumption  under  differ¬ 


ent  operating  conditions  for  the  gas  humidity.  This  should  also 

be  investigated  in  future  work. 
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